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ABSTRACT 

We carry out an analysis of a set of cosmological SPH hydrodynamical simulations of galaxy 
clusters and groups aimed at studying the total baryon budget in clusters, and how this 
budget is shared between the hot diffuse component and the stellar component. Using the 
TreePM+SPH GADGET- 3 code, we carried out one set of non-radiative simulations, and two 
sets of simulations including radiative cooling, star formation and feedback from supernovae 
(SN), one of which also accounting for the effect of feedback from active galactic nuclei 
(AGN). The analysis is carried out with the twofold aim of studying the implication of stellar 
and hot gas content on the relative role played by SN and AGN feedback, and to calibrate 
the cluster baryon fraction and its evolution as a cosmological tool. We find that both radia- 
tive simulation sets predict a trend of stellar mass fraction with cluster mass that tends to be 
weaker than the observed one. However this tension depends on the particular set of obser- 
vational data considered. Including the effect of AGN feedback alleviates this tension on the 
stellar mass and predicts values of the hot gas mass fraction and total baryon fraction to be in 
closer agreement with observational results. We further compute the ratio between the clus- 
ter baryon content and the cosmic baryon fraction, Yj,, as a function of cluster-centric radius 
and redshift. At i?5oo we find for massive clusters with A/500 > 2 x 10 14 h~ l M & that Yt is 
nearly independent of the physical processes included and characterized by a negligible red- 
shift evolution: I5.500 = 0.85 ± 0.05 with the error accounting for the intrinsic r.m.s. scatter 
within the set of simulated clusters. At smaller radii, i?25oo> the typical value of Yj, slightly 
decreases, by an amount that depends on the physics included in the simulations, while its 
scatter increases by about a factor of two. These results have interesting implications for the 
cosmological applications of the baryon fraction in clusters. 

Key words: cosmology: miscellaneous - methods: numerical - galaxies: cluster: general - 
X-ray: galaxies. 



1 INTRODUCTION 

Galaxy clusters, which stem from the collapse of density fluctu- 
ations involving comoving scales of tens of Mpc, are the largest 
gravitationally bound structures in the Universe. Since they are 
relatively well isolated systems, knowledge of their baryon con- 
tent is a key ingredient to understand the ph ysics of these ob- 
jects and their use as cosmological probes (see lAllen et alj|2~01 lb 
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iKravtsov and Borganill201 j . for recent reviews). The most massive 
galaxy clusters are of particular cosmological interest since their 
baryon content is expected to trace accurately the baryon content 
of the Universe. In the absence of dissipation, the ratio of baryonic- 
to-total mass in clusters should closely match the ratio of the cos- 
mological parameters measured from the cosmic micr owave back- 
ground (CMB) and therefore, M b /M to t ~ Ob/flm 1 White etalj 
1 19931 : lEvrardl I997l ; lEttori et all2 003; All en et al]2008l) . 



However, contrary to expectations, measurements of the 
baryon mass fraction in nearby clusters from optical and X- 
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ray observations have posed a challenge to this fundamental as- 
sump tion since they have reported smaller values than expected 
(e.g. lEttori et alj|2003l: iBiviano and SaluccfeOOrJ : iMcCarthv et"ai] 
12001) together with a possi ble intriguing trend with cluster mass 
(e.g., lLin etalj|2003l l2012h for a broad mass range of systems. 
With the precise measurement of the universal baryon fraction, 
ft = fib/Om = 0.167 ± 0. 004, from the Wilkin son Microwave 
Anisotropy Probe (WMAP-7: lKomatsu et aljboill) . these discrep- 
ancies have gained in physical importance and claim for dif- 
ferent explanations: physical processes which lower fh in clus- 
ters relative to the universal fraction (see e.g.. iBialek et al.|[200ll ; 
iHe et al.l 120061) . significant undetected baryo n components by 
standard X-ray and/ or optical techniques (see lEttori et al] 120031 ; 
iLin and Moh3 |2004). or a sy stematic underestimate of f2 m by 
WMAP ^McCarthy et al.ll2007h . 

In this sense, the correct determination of the gas mass fraction 
may be crucial. In fact, studies of the individual baryon components 
have shown that the stellar and gas mass fractions within RsotQ ex- 
hibit opposite behaviours as a function of the total system mass. 
In pa rticular, clusters have a higher gas mass f raction than groups 
(e.g.. lVikhlinin et alJl2006l:lArnaud et alj|2007l) , but a lower stellar 
mass fraction jLin et alj|2003l) . This has been interpreted as a dif- 
feren ce in the star formation efficiency between groups and clusters 
(e.g jDavid et al.ll 1990l : lLin et al.ll2003l : lLagana et alj2008l) . 

On the other hand, the mass dependence of the gas frac- 
tion and the discrepancy between the baryon mass fraction in 
groups/clusters and the WMAP value can be understood in terms 
of non-gravitational processes. In this regard, AGN heating, which 
can drive the gas outside the potential wells of host halos, can ex- 
plain the lack of gas within -R500 in groups. Therefore, groups also 
appear as critical systems to assess the universality of the baryon 
fraction and to understand complex physical processes affecting 
both the gas and the stellar components. 

The baryonic mass content of clusters consists of stars in clus- 
ter galaxies (satellite galaxies plus the brightest cluster galaxy or 
BCG), intra-cluster light (ICL, stars that are not bound to cluster 
galaxies), and the hot intra-cluster medium (ICM) whose mass ex- 
ceeds the mass of the former two stellar components by a factor 
of ~ 6. Therefore, to reliably estimate the cosmological param- 
eters from the baryonic-to-total mass ratio one should address all 
the baryonic components and not only the gas mass. In fact, to ex- 
plain the discrepancy between the observed baryon fraction and the 
cosmic value, the ICL is suggested to be one of the most important 
forms of missing baryons, acco unting for 6 — 22 per c ent of the total 
cluster light in the r-band (e.g., Gonzal ez et al.l20071) . or even more 
in cluster mergers (e.g.. |Pierini et alj|2008h . Cosmological simula- 
tions report that the ICL accounts for up to ~ 60 per cent of the to - 
tal of stars (e.g jMurante et al.ll2004l,l2007l ; |Puchwein et alj|2010h . 
However, although this large fraction of ICL is sufficient to explain 
the discrepancy with the cosmic value, it seems to be in contra- 
diction with observations. On one hand, this suggests that other 
mechanisms, like gas expulsion by AGN heating, may also be im- 
portant to account for the missing baryons. On the other hand, one 
should also bear in mind that the methods to identify ICL in sim- 
ulations and in observations are quite different, the former being 
often based on criteria of gravitational boundness of star particles, 



1 Ra ( A=2500, 500, 200) is the radius within which the mass density of a 
group/cluster is equal to A times the critical density (p c ) of the Universe. 
Correspondingly, Ma = A p c {z) (4 7r/3)iJ^ is the mass inside Ra- 



whil e the latter being ba sed on criteria of surface brightness limits 
(e.g. lRudick et al. II 20061) . 

The combination of robust measurements of the baryonic mass 
fraction in clusters from X-ray observations together with a deter- 
mination of fib from CMB data or big-bang nucleosynthesis cal- 
culations and a constraint on the Hubble constant, can therefore 
be used to measure f2 m (e.g., White and Frenkl 1991 ; White et al 



19931 : lEvrardl 19971: lAllen et al.ll2002l : lEttori et alj|2003l : ILin et al 



2003l : lAllen et al. 2003, 2004, 2008J). This method, remarkably sim- 



ple and robust in terms of its underlying assumptions, currently 
provides strong constraints on il m . On the other hand, measure- 
ments of the apparent redshift evolution of the cluster X-ray gas 
mass fraction, hereafter f g , can also be used to constr ain the ge- 
ometry of the Universe ( Sasaki 1996), its acceleration llPenll 19971) 



and, therefore, the dark energy equat i on of state (e.g., Ettori et all 
20031 12009|: I Allen et alj|2002l 120031 12004 iLaRoque et alj|2006l : 



Allen et alj|2008t 1201 ll) . This diagnostic exploits the dependence 



of the f g measurements (derived from the observed X-ray gas tem- 
perature and density profiles) on the assumed distances to the clus- 
ters, / g oc d , and relies on cluster baryon fractions being roughly 
universal and non-evolving over the redshift range where they can 
be observed (typically z < 1). Therefore, like Type la supernovae, 
massive clusters serve as standard calibration sources that test the 
expansion history of the universe. 

In this regard, early simulations by lEke et all dl998l) have 
been used to calibrate the depletion factor, i.e., the ratio by which 
the baryon fraction measured at ifesoo is depleted with respect 
to the universal mean (see, for instance, lAllen et al]|2008h . These 
non-radiative simulations indicate that, within the virial radius, the 
baryon fraction in clusters provides a measure of the universal 
mean that is only slightly biased low by ~ 10 per cent. The mag- 
nitude of this bias might change if additional physics is included in 
the simulations. Therefore, it is necessary and extremely useful to 
deepen in the analysis of how different implementations of bary- 
onic physics can affect the value of this bias, its evolution with red- 
shift, its radial dependence within clusters, or some statistics with 
massive galaxy clusters. 

Using non-radiative hydrodynamical simulations the expec- 
tation is that, for the largest (kT > 5keV), dynamically re- 
laxed clusters and for measurement radii beyond the innermost 
core (r > - R2500), fg should be approximate ly constant with 
redshift (e.g.. lEke et alii 19981 : Icrain et alj|2007l) . However, pos- 
sible systematic variations of f g with redshift can be accounted 
for in a straightforward manner if the allowed range of such vari- 
ations is constrained by numerical simulations or other comple- 



mentary data 1 Eke et alj 19981; Bialek et al. 


2001; Muanwong et al.l 


2002: Borgani et al.l2004l:lKav et al.l2004l: 


Ettori et al.l2004l2006l: 


Kravtsov et alj|2005: Nagai et al. 2007h. 



Therefore, it is clear that understanding the baryon mass frac- 
tion and its mass and redshift dependence is a crucial issue to 
understand better astrophysics in ga laxy clusters, e.g., the origin 
of the ICL (e.g. , IPierini et al.ll2008l) . star-formatio n history (e.g., 
Fritz et all 1201 ll) . metal-enrichment history (e.g., iKapferer et al .1 



20091) . the dynamical history of galaxy clusters, and the use of these 



system s to constrain the cosmological parameters (e.g. . lAllen et al .1 
l201lh . 

The purpose of the present work is to use a set of hydrody- 
namical simulations of galaxy clusters, characterized by different 
physical processes, to study how the fraction and spatial distribu- 
tion of baryons, as contributed both by the stellar component and 
by the hot X-ray emitting gas, are affected by the physical condi- 
tions within clusters and how these results compare with observa- 
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tions. We also analyse how the different baryonic depletions de- 
pend on redshift, baryonic physics, and cluster radius, determin- 
ing, therefore, some implications for the constraints on cosmo- 
logical parameters derived from gas mass fractions within clus- 
ters. In this regard, our analysis extends previous analyses of the 
baryon fraction in clust er simulations which included only non 



baryon traction in cluster simulations which included only non- 
radiative phy sics ( e.g., Evrardl 1990l : Metzler and Eyrardl 1994 



Navarro et alJll995J:lLubin et alJl996tlEke et al.lll998l ; lFrenk et all 
19991 : iMohr et al.ll 19991 ; iBialek et aljboOlh. the processes of cool- 



ing and star formation (e.g., iM uanwong et 



20041: 



Ettori etal 



Puchwein et al 
2012b . 



Sembolini et al.l 1201 



200 



ie processes ot cool- 
ail |2002|; iKav et al l 



20041 : iKravtsov et alj 120051 : lEttori et alj bool : 



, and t he effect of AGN feedback (e.g., 
l2010l : iFabian et alj I2OI0I : iBattaglia et all 



The paper is organized as follows: in Section 2, we describe 
our dataset of simulated galaxy clusters and the different physi- 
cal processes considered in re-simulating them; in Section 3, we 
present the results obtained from this set of simulations on the 
baryon, gas and stellar mass fractions as a function of cluster mass 
and we compare these results with different observational data sets; 
in Section 4 we calibrate the different baryonic depletions and anal- 
yse their dependences on redshift, baryonic physics and cluster ra- 
dius; and finally, in Section 5, we summarize and discuss our find- 
ings. 



2 THE SIMULATED CLUSTERS 
2.1 Initial conditions 

Our sample of simulated clusters and groups are obtained from 29 
Lagrangian regions, centred around as many massive halos iden- 
tified within a large-volume, low-res olution N-body cosmological 
simulation (see lBonafede et al]|201lL for details). The parent Dark 
Matter (DM) simulation followed 1024' 5 DM particles within a box 
having a comoving side of 1 h^ 1 Gpc, with h the Hubble con- 
stant in units of 100 km s _1 Mpc -1 . The cosmological model 
assumed is a flat ACDM one, with Q m = 0.24 for the matter 
density parameter, Qb = 0.04 for the contribution of baryons, 
Hq — 72 km s _1 Mpc -1 for the present-day Hubble constant, 
n s — 0.96 for the primordial spectral index and as — 0.8 for the 
normalisation of the power spectrum. Within each Lagrangian re- 
gion we increased the mass resolution and added the relevant high- 
frequency modes of the power s pectrum, following the zoomed ini- 
tial condition (ZIC) technique jTormen et al.lll997l) . Outside these 
regions, particles of mass increasing with distance from the tar- 
get halo are used, so as to keep a correct description of the large 
scale tidal field. Each high-resolution Lagrangian region is shaped 
in such a way that no low-resolution particle contaminates the cen- 
tral halo at z — at least out to 5 virial radifl As a result, each 
region is sufficiently large to contain more than one interesting halo 
with no contaminants within its virial radius. 

Initial conditions have been generated by adding a gas com- 
ponent only in the high-resolution region, by splitting each par- 
ticle into two, one representing DM and another representing the 
gas component, with a mass ratio such to reproduce the cos- 
mic baryon fraction. The mass of each DM particle is ?tidm = 



2 The virial radius, R v i r , is defined as the radius encompassing the over- 
density of virialization, as predicted by the spherical collapse model (e.g., 



lEkeet all 19961) . 



8.47 ■ 10 s /i _1 M and the initial mass of each gas particle is 
m Eas = 1.53 • 10 8 r'jtfp,. 



2.2 The simulation models 

All the simulations have been carried out with the TreePM- 
SPH GADGET- 3 code, a more e fficient version of the previous 
GADGET- 2 code dSpringelll2005t) . In the high-resolution region 
gravitational force is computed by adopting a Plummer-equivalent 
softening length of e = 511" 1 kpc in physical units below z = 2, 
while being kept fixed in comoving units at higher redshift. As for 
the hydrodynamic forces, we assume the minimum value attainable 
by the SPH smoothing length of the B-spline kernel to be half of 
the corresponding value of the gravitational softening length. 

Besides a set of non-radiative hydrodynamic simulations (NR 
hereafter), we carried out two sets of radiative simulations. 

A first set of radiative simulations includes star formation 
and the effect of feedback triggered by supernova (SN) explosions 
(CSF hereafter). Radiative cooling rat es are computed by fo llow- 
ing the same procedure presented by Iwiersma et al. I d2009h . We 
account for the presence of the cosmic microwave background 
(CMB) and of UV/X-ra y background radiation fro m quasars and 
galaxies, as computed bv lHaardt and Madavj J200 lh . The contribu- 
tions to cooling from each one of eleven elements (H, He, C, C, 
N, O, Ne, Mg, Si, S, Ca, Fe) have been pre-comp uted using the 
publi cly available CLOUDY photo-ionisation code l lFerland et alj 
1 19981) for an optically thin gas in (photo-)ionisation equilibrium. 
Gas particles above a given threshold density are treated as mul- 
tiphase, so as to provide a subresolution description of the inter- 
stellar medium, accord i ng to the model originally described by 
ISpringel and Hernquisll d2003al) . Within each multiphase gas par- 
ticle, a cold and a hot-phase coexist in pressure equilibrium, with 
the cold phase providing the reservoir of star formation. The pro- 
duction of heavy elements is described by accounting for the con- 
tributions from S N-II, SN-Ia and low an d intermediate mass stars, 
as described by iTornatore et al. I l l2007l) . Stars of different mass, 
distributed according to a Chabrier IMF dChabrieilPiool) . release 
metals o ver the time-scale determined by the mass-dependent life- 
times of IPadovani and Matteuccil dl993l) . Kinetic feedback con- 
tributed by SN-II is impleme nted according to the model by 
ISpringel and Hernquis tl fe003al) : a multi-phase star particle is as- 
signed a probability to be uploaded in galactic outflows, which is 
proportional to its star formation rate. In the CSF simulation set we 
assume v w = 500 kms -1 for the wind velocity, while assuming a 
mass-upload rate that is two times the value of the star formation 
rate of a given particle. 

Another set of radiative simulations is carried out by including 
the same physical processes as in the CSF case, with a lower wind 
velocity of v w = 350 kms" , but also including the effect of AGN 
feedback (AGN set, hereafter). In the model for AGN feedback, re- 
leased energy results from gas accretion onto supermassive black 
holes (BH). This model is based on th e original implementation of 
BH feedback bv lSpringel et alj 112005) (SMH), with some modifi- 
cations that will be described in detail by Dolag et al. (2012, in 
preparation). BHs are described as sink particles, which grow their 
mass by gas accretion and merging with other BHs. Gas accretion 
proceeds at a Bondi rate, and is limited by the Eddington rate. Once 
the accretion rate is computed for each BH particle, a stochastic 
criterion is used to select the surrounding gas particles to be ac- 
creted. Unlike in SMH, in which a selected gas particle contributes 
to accretion with all its mass, we included the possibility for a gas 
particle to accrete only with a slice of its mass, which corresponds 
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to 1/4 of its original mass. In this way, each gas particle can con- 
tribute with up to four "generations" of BH accretion events, thus 
providin g a more continuou s description of the accretion process 
(see also lFabian et al.ll2010i) . BH particles are initially seeded with 
a mass of 0.05 itidm.io, where m.DM,io is the DM particle mass in 
units of 10 10 ti -1 MQ. Seeding of BH particles takes place in ha- 
los when they first reach a minimum friend-of-friend (FoF) mass 
of 2.5 x 10 3 niDM.io (using a linking length of 0.16 in units of the 
mean interparticle separation in the high-resolution region), with 
the further condition that such halos should contain a minimum 
mass fraction in stars of 0.02. The latter condition guarantees that 
substantial star formation took place in such halos. This criterion 
prevents seeding BHs in halos possibly located at the border of the 
high resolution region, which spuriously contain a low amount of 
cooled gas, due to the interaction with nearby low-resolution DM 
particles. 

Eddington-limited Bondi accretion produces a radiated energy 
which corresponds to a fraction e r =0.1 of the rest-mass energy 
of the accreted gas, which is determined by the radiation efficiency 
parameter e r . The BH mass is correspondingly decreased by this 
amount. A fraction of this radiated energy is thermally coupled to 
the surrounding gas. We use e/ = 0.05 for this feedback efficiency, 
which increases to e/ = 0.2 when accretion enters in the quiescent 
"radio" mode and takes pla ce at a rate smaller than one-hundredt h 
of the Eddington limit (e.g. lSiiacki et all2007l ; lFabian et al.l2010l) . 



2.3 Identification of clusters 

The identification of clusters proceeds by running first a FoF al- 
gorithm in the high-resolution regions, which links DM particles 
using a linking length of 0.16 times the mean interparticle separa- 
tion. The center of each halo is then identified with the position of 
the DM particle, belonging to each FoF group, having the minimum 
value of the gravitational potential. 

Starting from this position, and for each considered redshift, 
a spherical overdensity algorithm is employed to find the radius 
Ra encompassing a mean density of A times the critical cosmic 
density at that redshift, Pc(z). In the present work, we consider 
values of the overdensit}0 A = 2500, 500 and 200. For the sake 
of completeness, we also consider the virial radius which defines a 
sphere enclosing the virial density A v i T (z)p c (z), predicted by the 
spherical collapse model (A v i r « 93 at z = and « 151 at z = 1 
for our cosmological model). In total, we end up with about 70 
clusters and groups having M V i r > 1 x 1O 14 /i _1 M0 at z = 0. 

Throughout this work all the quantities of interest will be eval- 
uated at the four different characteristic radii. Therefore, for each 
cluster, the hot gas, stellar, and baryonic mass fractions within a 
given radius Ra are defined, respectively, as 



/g(< 
/.(< 



M s (< R A ) 
M to t(< -Ra) 

M,(< Ra) 
M tot (< -Ra) 

M g (< Ra) + M,(< Ra) 
M t ot(< -Ra) 



(1) 
(2) 
(3) 



3 BARYON CONTENT OF CLUSTERS 

Figures |T| [2] and [4] show, respectively, the baryon, stellar, and gas 
mass fractions of our simulated clusters as a function of cluster 
mass M500- Only clusters with MsoosJ 3 x 10 13 h~ 1 M (3 within 
our NR, CSF, and AGN runs are shown. In order to compare with 
observational data, we u se some represen t ative observational sam - 

J2003h . iGonzalez et al. 



pies, mainly those f rom 



Lin et al 



iGiodini et al"] l l2009l) . rLagana et al 



fcOllh . and lZhang etaO 



2007) 



2011) 



3 The corresponding radii aproximately relate t o the virial radi us as 
(#2500, -Rsoo, -R200) ~ (0.2, 0.5, 0.7) il vir (e.g.. lEttori et alj2006h . 



For all the observational data sets we compare with, we compile 
in Table Q]fhe best fittings obtained for the baryon, gas, and stellar 
mass fractions. 

Before comparing our results with observational data, let us 
briefly describe the main properties of the different observational 
samples (for further details, we refer to their corresponding papers). 
Knowing how the observational data have been derived is important 
to understand not only discrepancies between our simulations and 
the observations but also the differences between the observational 
results. 

lLagana et all d201 lh and IZhang et all d201 lh investigate the 
baryon mass content for a subsample of 19 clusters of galaxies ex- 
tracted from the X-ray flux-limited sample HIFLUGCS. For these 
clusters, the above authors measure total masses and characteristic 
radii on the basis of a rich optical spectroscopic data set, the physi- 
cal properties of the intra-cluster medium using XMM-Newton and 
ROSAT X-ray data, and total (galaxy) stellar masses utilizing the 
DR-7 SDSS multi-band i maging. Using gas m ass measurements 
from X-ray observations, lLagana et al.l d201 lh use a scaling rela- 
tion between the gas and the total mass t o determine t he tot al clus- 
ter mass. Following a different approach. lZhang et al.l d201 lh derive 
cluster masses from mea surements of the "har monic" velocity dis- 
persion as described by iBiviano et al .1 d2006l) . In both studies, to 
obtain the contribution of galaxies to the stellar mass, they use the 
optical data for selected member galaxies within -R500 to compute 
their luminosity function in the i-band, performing an analytical fit 
using two Schechter functions. Then, they adopt dif ferent mass-to- 
light r atios for ell ipticals and spir als, taken from Kau ffmann et al .1 
d2003h assuming a lSalpeteil ( ll955f) IMF, to compu te the stellar mass 
from the optical luminosity. In Lin et al .1 d2003l) they use the ob- 
served X-ray mass-temperature relation together with published 
X-ray emission weighted mean temperatures, 2MASS second in- 
cremental release NIR data, and X-ray imaging data to explore 
trends in the NIR and X-ray properties of a sample of 27 nearby 
galaxy clusters. The total mass of the clusters in this sample is es- 
timated from an observed M500 — Tx relation. The stellar masses 
are obtained from the cluster luminosity function as derived from 
the magnitudes in the K s band. Using a Schechter function they 
estimate the total cluster luminosity and, finally, they obtain the 
total stellar masses, as contributed by satellite galaxies and BCG, 
by multiplying the total luminosity by the average stellar mass-to- 
light ratio for each cluster which takes into account the varying 
spiral galaxy fractio n as a function of the cluster temperature. In 
IGiodini et all d2009l) . 9 1 candidate X-ray groups/poor clusters at 
redshift 0.1 ^ z < 1 are selected from the COSMOS 2 deg 2 sur- 
vey, based only on their X-ray luminosity and extent. They use X- 
ray detection, gravitational lensing signal, optical photometric and 
spectroscopic data of the clusters and groups identified in the COS- 
MOS survey. Total cluster masses are derived from a Lx — M500 
relation. The stellar mass of a galaxy is obtained from the conver- 
sion of the AVband luminosity using an evolving galaxy-type de- 
pendent stellar mass-to- K s -band luminosity ratio. Since this sam- 
ple is mostly composed of groups, it is complemented with the 27 



Baryon Census in Hydrodynamical Simulations of Galaxy Clusters 5 



nearby X-ray selected clusters in the sample bv lLin etalj J2003h . 
where the total and stellar masses are derived in a consistent man- 
ner. To reduce sy stematic effects, they use the scaling relations in 
IPratt et all d2009l) . based on hydrostatic mass estimates, to derive 
the total gas fractions in both samples. The total sample of 118 
groups and clusters with z ^ 1 spans a range in M5 00 of ~ 10 13 - 
10 15 M Q . On the other hand. lGonzalez et alj J2007h use an optical 
sample of 24 nearby clusters and groups for which they obtain drift 
scan imaging in Gunn i using the Great Circle Camera on Las Cam- 
panas 1 m Swope telescope. This sample is composed of systems at 
0.03 ^ z sj 0.13 that contain a dominant BCG. To obtain the total 
masses and cluster radii they derive calibrations o f the a — R500 
and a — M500 relations using the clusters from IVikhlinin et al.1 
d2006l) that also have published velocity dispersions. They deter- 
mine the luminosity of the BCG+ICL component by fitting the 
surface brightness distribution in each cluster out to a radius of 300 
kpc from the BCG. Then, they use the separate r 1//4 best-fit profiles 
of these two components to build a two-dimensional model image 
from which they determine the flux within a given circular aper- 
ture. On the other hand, the luminosity of the cluster galaxies lying 
within the same aperture is computed by summing the flux of all 
galaxies fainter than the BCG and brighter that mi = 18. Then, 
they use a relation for the luminosity dependence of the mass-to- 
light ratio in the I band to convert from total luminosity to total stel- 
lar mass. Since they lack measurements of the mass of hot gas in the 
ICM for this sample, they fit the behaviour of the stellar mass frac- 
tion with cluster mass and use this relation to derive total baryon 
fractions for clusters with published X-ray gas fractions. We also 
compare our results for the gas mass fractions with a sample of ob- 
served groups and clusters at z ^ 0.2 selected from the X-ray sarn- 
ies o f lvikhlinin et alj d2006h . lArnaud et al] d2007h and lSunetaf] 
20091) . These authors computed gas mass fractions at R500 from 
hydrostatic mass estimates for a combined sample containing 41 
systems within a range of masses of [1.5x 10 13 ,1.1 x 10 15 ] Mq. 

Due to the different observational methods used to derive the 
main cluster properties, some differences are expected between the 
baryon census provided by these samples. In addition, it is neces- 
sary to point out tha t , whe n comparing the stellar mass fractions, 
only iGonzalez et al.1 d2007h take into account the controbution of 
the ICL component. In our case we also consider the total (galax- 
ies+ICL) stellar contribution within clusters. In the following, after 
comparing simulation results to the observed total baryon budget 
in clusters and groups, we will dissect the separate contribution of 
stars and of hot X-ray emitting gas. 



3.1 Baryon mass fraction 

As shown in Fig.Q] in our NR simulations the baryon mass fractions 
within -R500 is nearly independent of cluster mass and is systemat- 
ically lower than the assumed cosmic value by J; 10 per cent. This 
result is cons istent with previous analyses, also based on SPH sim- 
ulations (e.g. lEke et al ] |l998l ; lEttorietalJl2006h . which also found 
a comparable undere stimate in the cluster b aryon fraction. Using a 
Eulerian AMR code, iKravtsov et alj J2005h also measured a clus- 
ter baryon fraction in non-radiative simulations below the cosmic 
value, although in their case the underestimate was of about 5 per 
cent within -R500 • A similar behaviour is also found for our radia- 
tive CSF simulations, thus indicating that the processes of star for- 
mation and galactic winds triggered by SN explosions determine 
the fraction of baryons to be converted into stars, without however 
changing the overall baryon budget within i?soo- 

As we include the effect of BH feedback in the AGN simula- 
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Figure 1. Baryon mass fraction as a function of cluster mass M500 ■ Results 
from our NR, CSF, and AGN runs are represented by black circles , blue tri- 
angles , and red stars, respectively . The observational samples from lLin et alj 
J2003I) . and lLagana et al. (201 1) are shown as shaded regions in gr een and 
orange, respectively, whereas the sample from lGiodini et al] 120091) is rep- 
resented by a grey-striped area. These regions correspond to the best fits 
obtained together with their corresponding errors (see Table [T}. The hori- 
zontal continuos line stands for the cosmic baryon fraction assumed in our 
simulations. 



tions, there is a significant baryon depletion in poor clusters and 
groups, whereas results are nearly the same as for the NR and CSF 
cases for Af 5 oo~ 2 x 10 14 h^ 1 Mq. This result of a decreasing 
baryon fract i on at lo w masses is in l ine wi th tho se presented by 
iFabian et all J2010h . IPuchwein etatl J2010h . and iMcCarthv et~al] 
J201 1), who also included the effect of BH feedback in their simu- 
lations of galaxy groups and clusters. This effect of baryon deple- 
tion within groups witnesses the efficiency that BH feedback has 
in displacing gas outside forming halos. This effect takes mostly 
place at relatively high redshift, z ~ 2 — 3, around the peak of 
the BH accretion efficiency. At these epochs, the energy extracted 
from BHs increases the gas entropy to levels such to prevent this 
gas from being subsequent ly re-accreted within group-sized halos 
(e.g. lMcCarthv et alj20llh . 

As for the comparison of simulation results with observations, 
it is quite remarkable that a good agreement is only achieved for the 
AGN model. This result confirms that a feedback mechanisms only 
based on SN explosions can not be responsible for the decreasing 
trend of the baryon budget within halos of decreasing mass. 



3.2 Stellar mass fraction 

We show in Fig. [2] the stellar mass fraction in our radiative runs. 
As expected, the effect of including AGN feedback is that of re- 
ducing the stellar content of galaxy systems by about 30 per cent, 
nearly independent of cluster mass. As for the comparison with 
previous simulation results, we note that the clusters simulated by 
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Figure 2. Stellar mass fraction as a function of cluster mass M500 ■ Results 
from our radiative simulations, CSF and AGN, are represented b y blue tri- 
angles and red stars, respectiv ely. The observational sa mples from lLin et alj 
120031) . iGonzalez et al.l 120071) and lLagana et alj 1201 ll) are shown as shaded 
regions in green, cyan and orange, respectively. These regions represent the 
best fits obtained from each observational sample as compiled in Table [T] 
The horizontal continuos line stands for the assumed baryon mass fraction 
in our simulations. 



iPuchwein et al with AGN feedback have stellar fractions, 

which are larger by about a factor of 2 than the stellar fractions 
found in our AGN simulation. We can understand this difference by 
keeping in mind that the amount of stars formed in simulations de- 
pends rather sensitively on how the SN feedback i s included (e.g., 
ISpringel and Hernquislll2003btlBorg ani et al.ll2006r) . In the simula- 
tions bv lPuchwein et alj J2010t) . all the feedback from star forma- 
tion was injected thermally, without including kinetic SN feedback 
as we did in our simulations. As these authors also noticed (see their 
Sect. 3.1), using kinetic feedback, in addition to thermal feedback, 
can significantly reduce the amount of stars formed in their simu- 
lations by a factor of 2, resulting, therefore, in a good consistency 
with the stellar mass fractions obtained in our AGN runs. 

As for the comparison with observational results, we find that 
our CSF simulations produce a too large stellar fraction in mas- 
sive galaxy clusters, independent of the observational data set we 
compare to. While simulations with AGN feedback give results 
closer to observations, the level of agreement is quite sensitive to 
the observational r esult we refer t o. For instance, a comparison 
with the results bv lGonzalez et alj J2007I) would imply that in no 
case simulations reproduce the steep mass dependence of / * , inde- 
pende ntly of the feedback mechanism included (see also lAndreorJ 
l2O10h . On the other hand, a closer agreement with observations 
w ould be obtained fro m Fig. [2] by referring instead to the results 
bv lLagana et al J 1 20 1 1 ). The inclu sion of the ICL component in the 
analysis bv lGonzalez et alj d2007t> could explain part of the differ- 
ence with respe ct tolLagana et alj j2oTTb . although apparently not 
all of it (see also lZhang et alj201 lh . Clearly, some caution must be 
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Figure 3. Fraction of the stellar mass found in the BCG+ICL component as 
a function of cluster mass M500 ■ The upper panel shows the values obtained 
from our radiative runs without AGN feedback, the CSF runs, while the 
lower panel displays the results obtained from our runs including AGN. We 
compare our re s ults w ith the observed BCG+ICL luminosity fractions from 
lGonzalezet"a?]fe007t) . 



used when comparing observational and simulated samples, ow- 
ing also to the different approaches followed by different authors 
to measure stellar mass fraction from data, the depe ndence of the 
inferred stellar mass on the choice of the IMF (e.g.. lLagana et al] 
I201 it iLeauthaud et alj|2012h . and systematic uncertainties in the 
measurement of the total cluster mass. 

In order to better understand how much intra-cluster stars con- 
tribute to the total stellar mass budget in our simulations, it is im- 
portant to distinguish between the stellar content of the BCGs, of 
the satellite galaxies and of the ICL component. Aside from some 
exceptions, the central galaxy in a cluster, which is the closest to the 
minimum of the cluster potential well, is typically also the brightest 
cluster galaxy. Therefore, for simplicity, we will use the abbrevia- 
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tion BCG when referring to the central galaxy of a simulated clus- 
ter or group. Due to its low surface brightness, observations of the 
ICL component, which is a smoothly distributed stellar component 
typically peaked around the BCG but extended to larger radii, are 
difficult, resulting in a significant uncertainty in the current obser- 
vational constrain t s of the amount of IC L present in clusters (e.g., 
IZibetti et al]2005t IGonzalez et al.ll2007l) . 

While identifying member galaxies within galaxy clusters in 
hydr odynamical simulat ions is a relatively straightforward task 
(e.g., lOnions et al.ll2012h . distinguishing between stars in the dif- 
fuse stellar component and in the central galaxy is not so trivial. In 
order to do so, we use a modified version of the SUBFIND algo- 
rithm ( Springel et al. ll200ll ; lDolag et alj|2009l) . In its original ver- 
sion jSpringel et al .11200 lh . SUBFIND identifies star particles that 
are associated with satellite galaxies residing within a cluster-sized 
halo. All the star particles not associated to satellite galaxies are as- 
signed to the central galaxy of the main halo, without distinguishing 
between those associate d to the actual BC G and those belonging 
to the surrounding ICL. iDolag et all j2009h pointed out that BCG 
and ICL stars show different phase-space distributions and imple- 
mented this property in SUBFIND, making it able to distinguish 
among the different stellar components. For more details about this 
modified version o f SUBFIND we refer the reader to the work by 
lDolagetal.|j20ld) . 

An intrinsic difficulty in properly comparing observations and 
simulation results on the amount of ICL is due to the intrinsically 
different procedures usually adopted to identify diffuse stars in real 
and in simulated data. While observations generally use a criterion 
based on surface brightness limit, ICL analysis in simulations is 
generally based on identifying star particles that are not gravitation- 
ally bound to galaxies. In addition, very faint ICL component can 
not be detected in observations while it is present in simulations. 
While we defer to a future paper a homogeneous comparison be- 
tween intra-cluster light in observations an d simulations, we show 
here a comparison between the results by iGonzalez et al. I j2007l) 
on the amount of stars present in BCG and ICL, and correspond- 
ing simulation results. In fact, considering the total stellar content 
of BCG and ICL overcomes at least the ambiguity in the surface 
brightness limit below which the BCG halo has to be considered 
as part of the BCG. This also avoids choosing among th e different 
ICL definitions in the literature (e.g. IZibetti et alj|2005r) allowing, 
therefore, a more straightforward comparison between our results 
and other simulations and observational studies. In Fig. [3] we plot 
the fractions of stellar mass found in the BCG+ICL components 
in our simulated sample of clusters with M500 ^ 1-5 x 10 Mq. 
Results obtained from our radiative simulations without and with 
AGN feedback (CSF and AGN runs), are shown in the top and bot- 
tom panels, respectively. We compare our data wi th the observa- 
tional constraints on the BCG+ICL fraction from IGonzalez et all 
j2007l) , shown as blue triangles with error bars. 

Although our results in general confirm a decreasing trend 
with cluster mass of the fraction of stars contributed by BCG and 
ICL, they predict a too large v alue of this frac t ion in comparison 
with the observational result bv lGonzalez et alj j2007l) . In the CSF 
simulations we obtain BCG+ICL fractions of roughly ~ 60 per 
cent for massive clusters, and of about ~ 90 per cent for groups. 
These values are larger than those observed by IGonzalez et al .1 
j2007h . especially for massive clusters. However, we have to take 
into account that there is some uncertainty in the mass assigned to 
the BCG and ICL components due to the analysis metho d used for 
separating them (see, for instance. IPuchwein et al.ll2~010l) . In addi- 
tion, given that the distribution of the BCG+ICL component is more 
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Figure 4. Gas mass fraction as a function of cluster mass M500 . Results 
from our NR, CSF, and AGN runs are represented by black circles, blue tri- 
angles and red stars, respectively. We compare our results with two differ- 
ent observational samples: a combined sam ple of 41 clusters and groups 



from IVikhlinin et all 12006b . lArnaud et all j2007t) and ISun et alj j2009l 
(V06+APP07+S09), shown as the ora nge region, and the s ampl e obtained 
from the combination of the data by Izhang et alj 1201 ll) and ISun et all 
42009b (Z11+S09), shown as the green area (see Table [TV The horizontal 
continuos line marks the cosmic value of the baryon mass fraction assumed 
in our simulations. 



concentrated than the distribution of the satellite galaxies, these val- 
ues are sensitive to the radius inside which they are measured. 

When AGN feedback is included, we find an even larger frac- 
tion of stars in the BC G+ICL component, a re sult that is consistent 
with that presented bv lPuchwein et al] j201fj|) . The reason for this 
result is that, although the stellar mass of the BCG+ICL component 
decreases when including AGN feedback, the total stellar mass de- 
creases even more (see Fig.|2). A relative increase of stars in BGC 
and ICL is the consequence of the combination of two different ef- 
fects. On the one hand, the effect of AGN feedback is mainly that of 
truncating the star formation of clusters a t high redshift, zj; 3 (e.g. 
iFabian et al.ll2010l : lMcCarthv et alizOl ]]). On the other hand, most 
of the dynamical origin of the ICL is a ssociated with the assembly 
of the BGC (e.g. iMurante et ai]|2007l) . Since mergers continue to 
take place after star formation is quenched by AGN feedback, they 
keep unbounding stars from galaxies into the diffuse intra-cluster 
components. Since this process is not compensated by fresh star 
formation in the presence of AGN feedback, the net effect is that of 
increasing the fraction of stars that end up in the ICL. In a future 
analysis (Cui et al., in preparation) we will carry out a more de- 
tailed comparison of ICL inventory and properties in simulations, 
by reproducing in their analysis the same criteria to identify ICL as 
in observational data. 
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3.3 Gas mass fraction 

As shown in Fig. [4] including radiative physics in the simulations 
has the e xpected effect of decreasing the gas fraction within -R500 
(see also iKravtsov et al .1 120051: iFabian et alJIgOlCl I Puchwein et al.l 
l2010l : lMcCarthv et al.ll201 lUSembolini et alj|2012l) . by an amount 
which is more pronounced in poor clusters and groups. As for the 
effect of including different feedback mechanisms, a comparison 
between our simulations with and without AGN feedback shows 
that the two feedback schemes predict rather similar values of the 
gas fraction at the mass scale of groups while simulations includ- 
ing AGN feedback predicts slightly more gas within large clusters. 
Clearly, the similar values of f g in groups do not imply that feed- 
back does not have any effect on such systems. In fact, a compar- 
ison with Figs. Q] and [2] highlights that AGN feedback tends to re- 
move baryons from the potential wells of galaxy groups. At the 
same time, suppression of star formation partially prevents removal 
of gas from the hot diffuse phase within R500, thereby acting as a 
compensating effect such that the resulting gas fraction turns out 
to be similar for the two feedback schemes. As for higher-mass 
halos, AGN feedback becomes less efficient in removing baryons 
from halos (see also Fig. [T}, so that suppression of star formation 
causes a slightly larger fraction of baryons to remain in the diffuse 
phase, so that f g in this case increases as a result of a more effi- 
cient feedback. This differential effect of AGN feedback in low- 
and high-mass halos is generally quite weak, although it goes in 
the direction of better reproducing the observed trend of f g with 
halo mass. 

From the analysis of Fig. [4] we conclude that, in general, our 
results on the values of f g at the scale of rich and poor clusters, es- 
pecially in the presence of A GN feedback, ar e in line with the ob- 
servat i onal results obtain ed b ylVikhlinin et al . ( 2006), Arnaud et al .1 
d2007hJSun et~aIU2009l) . and lZhang et all d20 1 lh ■ 



4 CALIBRATION OF THE BARYONIC BIAS 

After having compared simulation results on the different baryonic 
components with observational data, in this Section we use our re- 
sults to calibrate the different baryonic depletions and to analyse 
their dependences on redshift, baryonic physics and cluster radius. 

For the sake of comparison with previous works, we de- 
fine the gas, stellar and baryon depletion factors (from now on 
Yg, Y*, and Yb, respectively) as the ratios between f s , /„ and 
/b = f g + /*, and the cosmic value adopted in the present 
simulations, fib/fi m = 0.167. Accordingly, we should measure 
Yb = 1 within clusters as long as they are fair containers of cos- 
mic baryons. Any deviation from this value has to be interpreted 
as due to the presence of a "baryonic bias", whose origin can be 
due either to gas dynamical effects at play during the hierarchi- 
cal assembly of clusters, or to star formation and feedback effects 
that causes sinking or expulsion, respectively, of baryons from the 
cluster potential wells. Th e non-radia t ive sim ulations of hot, mas- 
sive c lusters published by lEke et al" I J 19981) (see also ICrain et all 
|2007|) give Y M = 0.83 ± 0.04 at #2500, and are consistent with 
no redshift evolution of Yb for z < 1. Neverthe less, simulations 
including different models of baryonic physics JKay et al. I l2004l ; 
lEttori et alj|2006l : ICrain et al.ll2007l ; lNagai et alj|2007l) allow for a 
range of evolutions. We note, however, that these previous analyses 
either lack sufficient statistics of massive systems, which are rele- 
vant for cosmological applications, or the inclusion of an efficient 
feedback mechanism, like that provided by AGN, which provides 



a realistic description of star formation in the central regions of 
galaxy clusters. 

The results that we will present in the following are relevant 
to test the robustness of the calibration through simulations of the 
baryon bias, i.e. of the deviation of the baryon content of clusters 
from the cosmic value, that one needs to correct for in the cosmo- 
logical application of the gas mass fraction. 

4.1 Radial dependence of the baryonic bias 

We show in Fig. [5] the mass dependence of Y g and Yb at 
2 = (up and bottom row, respectively) within the two 
characteristic radii i?soo (left column) and -R2500 (right col- 
umn). As for -R500, it typically corresponds to the most exter- 
nal radius out to which detailed X-ray observations, possibly 
in combination with Sunyaev-Z eldovich (SZ) observations (e.g. 
IPlanck Collaboration et al.l201ll) , allow one to trace the gas content 
within clusters, while R2500 is the typical radius within which gas 
content is traced for distant clusters, when using the evolution of 
the ga s fraction in clusters as a cosmological probe (e.g. l Allen et al .1 
l2008h . Therefore, for these two radii, we show the mean values of 
the gas and baryon depletion factors obtained for our sample of 
simulated clusters binned in five linearly equi-spaced bins in M500 • 
All clusters with M500K1 1 X 10 13 /i _1 M have been considered. 
The mean values within each mass bin are shown along with error 
bars representing one standard deviation within the corresponding 
mass interval. 

The left panel of this figure summarizes the simulations re- 
sults shown in Figs.[T]and|4] Using the mass binning, it is now more 
clear that the depletion in baryon content within R500 is more pro- 
nounced and with a stronger mass dependence for the simulations 
including AGN feedback, at least for low-mass systems. As already 
discussed, the larger baryon depletion in clusters simulated with 
AGN feedback is the result of the efficiency of this feedback mech- 
anism in removing baryons from the potential wells of forming 
groups at redshift z ~ 2-3, around the peak of gas accretion onto 
SMBHs. On the other hand, for masses M 500 > 2 x 10 14 h~ 1 M Q 
we find that the baryon fraction within R500 underestimates the 
cosmic value by about 15 per cent, nearly independent of mass 
and of the physics included in the simulations. The r.m.s. disper- 
sion around this values is of about 3 per cent for the NR and CSF 
simulations, which increases to about 5 per cent for the AGN sim- 
ulations. This result for Yj, is different from the behaviour of gas 
depletion, which shows no flattening for high-mass systems. Fur- 
thermore, values of Y g for the AGN simulations are systematically 
larger than for the radiative simulations including only SN feed- 
back, as a result of the suppressed star formation in the former case. 

These results suggest that a mass-independent correction can 
be calibrated from simulations to infer the cosmic baryon fraction 
from the corresponding quantity derived for massive clusters, a cor- 
rection that is likely independent of the uncertain knowledge of the 
physical processes at play in the ICM. However, these results also 
highlight that accurately recovering the baryon fraction from gas 
mass measurements involves accurately accounting for a correc- 
tion associated to stellar mass, which generally depends on cluster 
total mass. 

Results on gas and baryon depletions are somewhat different 
within J?2500 (right column of Fig. |5). In this case, both Y g and 
Yb show as steady increase with cluster mass, with no evidence 
for a flattening at high masses, and a larger intrinsic scatter in 
their values. Furthermore, a small but sizeable dependence of Y, 
on the physics included in the simulations exists even for the high- 
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Figure 5. Mass dependence of the gas depletion factor Y g (upper row), and total baryon depletion factor Yb (lower row), computed at -R500 and -R2500 (left 
and right columns, respectively). Results are shown for the set of simulated clusters with Afsoo~ 1 X 10 13 /i -1 Af0 identified at z = within our NR, CSF, 
and AGN runs. Clusters are binned in five linearly spaced mass bins. Different line types represent the mean values obtained within each mass bin for each of 
our simulations, while error bars stand for Icr standard deviations computed within the subset of clusters within each mass bin. For reasons of clarity, lines 
corresponding to the different physical models have been slightly displaced along the x-axis. 



est mass systems. Quite interestingly, the largest values of Yj, are 
obtained for the CSF simulations, as a result of the strong over- 
cooling, not efficiently counteracted by the SN feedback, which 
causes a large amount of baryons to condense in the central halo 
regions. On the other hand, the smallest Yj value is obtained in the 
presence of AGN feedback, which is quite efficient in displacing 
gas from central regions even for the most massive clusters. In gen- 
eral, these results indicate that the baryon depletion within -R2500 
is more sensitive to the detailed description of the feedback process 
which regulates the cooling-heating cycle. As a result, care must 
be taken in the use of simulations to exactly calibrate the correc- 
tion for baryon depletion to observations providing gas and baryon 
fractions at such smaller cluster-centric radii. 

As shown in Fig. [5] the population of the most massive clus- 
ters is characterized, especially within Rsoo, both by a remark- 
ably small intrinsic scatter in their baryon budget, and by a stabil- 
ity against the different physical descriptions of the ICM. As such 
these massive clusters are those which can be more reliably cali- 
brated for cosmological applications of the cosmic baryon fraction 
test. Therefore, in the following we will restrict the analysis of these 
radial distributions only to clusters with M500 > 2 x 10 14 h _1 MQ. 
Within our simulations we identify about 40 of such objects at 
z — 0, a number that reduces to 10 at z = 1. 



The effect of the different physical models considered in our 
simulations on the distribution of baryons can be better understood 
from the analysis of the radial distribution of the different baryonic 
components within clusters. Figure[6]shows the mean radial distri- 
bution out to 4 -R500 of the baryonic, gas and stellar depletions at 
z — (left column) and z = 1 (right column) for our subsample 
of massive clusters within each of the physical schemes adopted in 
our simulations. 

Regardless of the baryonic processes included in our re- 
simulations, the baryonic depletion at z = for radii r/i?soo ^ 
0.4 approaches a value of ~ 85 per cent of the cosmic value, 
showing similar values at z = 1 but with larger scatter. This 
baryonic depletion starts to converge to the expected value at 
r ~ 3 x . Rsoo, consistent with results found in previous simu- 
lations (e.g jEke etalll 19981 ; lEttori et alj|200d) . In these outer re- 
gions of clusters, the gas mass dominates the baryon budget. In 
general, the gas depletion increases from inner to outer regions 
and shows slightly higher values at low redshifts. On the contrary, 
the stellar depletion decreases when moving towards more exter- 
nal regions and, therefore, if we move towards more internal radii 
(r/Rsoo ^ 0.1) the stellar mass clearly dominates the baryon con- 
tent in the radiative runs. In these central regions, the non-radiative 
simulations produce lower values of the baryonic depletion than 
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Figure 6. The mean radial profiles out to 4 X -R500 of the baryonic, gas and stellar depletions (from top to bottom panels) at 2 = (left column) and z = 1 
(right column) for the subsample of massive clusters with Afsoo > 2 X 10 14 h~ 1 Mq. In each panel, mean profiles for the NR, CSF and AGN simulated 
clusters are shown with the continuous black, dashed-dot blue, and long-dashed red lines, respectively. For each model, dotted lines around the mean profiles 
indicate the region corresponding to 1<t standard deviation around the mean. Within each panel and from left to right, dotted vertical lines indicate the position 
of the mean value of i?2500, R200 and R v i r (in units of Rtmo ) for the sample of clusters for which the profiles are computed. 



the radiative runs which are, indeed, characterised by a steep in- 
ner slope. When cooling and star formation are included, the gas 
can cool and form stars and, therefore, it sinks deep into the poten- 
tial wells of the clusters. If AGN feedback is also added, its main 
effect is that of heating the surrounding gas producing, therefore, 



smaller baryon mass fractions in the cluster center. In general, the 
different baryonic depletions obtained in the inner regions of clus- 
ters from the CSF and AGN simulations are comparable with each 
other. However, whereas at 2 = 1 the baryon and stellar depletions 
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are comparable or even slightly higher in the AGN runs, the roles 
are clearly inverted at z = as a consequence of the AGN heating. 



4.2 Redshift evolution of the baryonic bias 

Besides assessing the stability of the baryon and gas depletions at 
z — 0, measuring the evolution of such quantities is also required 
for measurements of the gas fraction over a large redshift baseline 
to be used to recover the redshift-distance relation. A comparison 
of the profiles of baryon and gas depletion computed at z — 
and z — 1 (see Fig. |6) shows that outside the central regions this 
evolution is generally rather mild. 

To quantify this evolution, we compute the values of the de- 
pletion factors at different redshifts, z = 0, 0.3, 0.7, 0.8 and 1. 
At each redshift, the analysis is done only on those clusters having 
mass M 5 oo > 2 x 10 14 h~ 1 M & . We show the redshift evolution 
of these quantities at -R500 and -R2500 in the left and right panels 
of Fig. [7] respectively, along with their respective intrinsic scatters, 
given by the error bars. A compilation of these values, for the dif- 
ferent redshifts and simulation sets, are also reported in Table[2] 

In order to parameterize a possible evolution of the values of 
the depletion factors, we use the expression 

Yi(z) = Y ,i(l + a Yi z), (4) 

where the subscript i is equal to 6 or g when referring to the total 
baryon or gas content, respectively. The values of the parameters 
Yb,i and ay ; are computed through a \ 2 minimization procedure, 
with the weights of the data points reported in Table [2] provided by 
their corresponding intrinsic scatter. 

In Table[3]we report the best-fitting values of these parameters 
obtained for each simulation set, within different radii of interest, 
from i?25O0 out to the virial radius R v i r . 

The results displayed in Fig. [7] show that, independently of 
the considered radius or physics, the baryon depletion factor does 
not evolve significantly with redshift, at least since 2 = 1. Within 
^2500 (right panel) the dissipative action of radiative cooling in the 
CSF runs slightly increases the average value of Vj, with respect to 
the AGN simulations, bringing it very close or even above to that 
of the non-radiative simulations, with Yb ~ 0.85, constant across 
the considered redshift range. On the other hand, the presence of 
AGN feedback is effective in preventing gas from accreting onto 
the central regions, thus decreasing the baryons fraction to Yj ~ 
0.80, also independent of redshift. 

As for results at -R500, we find a smaller scatter and much 
better agreement among the different physical models, thus high- 
lighting that the different physical descriptions of the ICM have 
a negligible impact on the total amount of baryons at such larger 
cluster-centric radii. At such radii, we find Yb — 0.85 virtually in- 
dependent of redshift, with some departure for the AGN simulations 
at z = 1, probably due to the limited statistics of massive clusters 
at the highest considered redshift. Therefore, a sizeable decrease 
in the baryon fraction when moving inwards to -R2500 is detected 
when including the more realistic feedback scheme based on the 
effect of AGN. 

As for the gas mass fraction, the inclusion of radiative physics 
decreases its value with respect to the non-radiative simulations, 
both at -R2500 and at i?soo- As expected, this decrease is more 
pronounced at smaller radii and for the simulations only includ- 
ing the effect of SN feedback. As for the AGN simulations, we find 
Y g ~ 0.5-0.6 within R2500, quite independent of redshift, with a 
significant scatter, ay g — 0.1, over the whole range of redshift. 



This value increases to Y g ~ 0.6-0.7 within R500, also nearly con- 
stant in redshift, but with a reduced intrinsic scatter of <ry g ~ 0.05. 
The behaviour obtained for the CSF simulations is pretty similar 
but with lower values for the gas fraction: within -R2500, Y g ~ 0.5, 
whereas it is Y g ~ 0.6-0.7 at R500. Quite remarkably, in all cases 
such values are nearly independent of redshift. 

In general, our results for the NR case are i n agreement with 
tho se from non-radiati ve simulations presented bv lEke et alj i ll 9981) 
and lEttori et alj ([2006), both based on SPH simulations, while they 
are slightly, but s ystematically, lower b y about 5 per cent than 
those obtained by iRravtsov et alj J2005T) from AMR simulations. 
Although this difference is quite small, it is still comparable to, 
or larger than, the difference induced by the presence of different 
physical processes in simulations. Although it remains to be seen 
whether such a difference between predictions of SPH and AMR 
codes persists when including radiative physics, its presence warns 
on the need of understanding in detail the performances of different 
hydrodynamical methods in the calibration of the gas mass fraction 
test through simulations. In general, our results on non-radiative 
simulations including only SN feedback are in line with those pre- 
sented by other authors (e.g. Muanwong et al. 2002; Kr avtsov et alj 
l2005l : lEttori et alj2006l : ISembolini et alj2012l) . However, while the 
comparison between results from different non-radiative simula- 
tions is relatively straightforward, when extra-physics is included 
the results on the distribution of the different baryonic components 
are sensitive not only to the nature of the feedback sources included 
(i.e. SNe vs. AGN), but also to the details of the numerical imple- 
mentation (i.e. thermal vs. kinetic feedback, dependence of cooling 
rates on local gas metallicity, numerical resolution). These aspects 
have to be taken into account when performing such a comparison 
between results from different authors. 

In principle, the results of our analysis can be used to set 
priors on the parameters which determine the amount of gas de- 
pletion within a given aperture radius and its redshift evolution, 
when deriving cosmological parameters from observations of the 
baryon fraction in massive clusters. Overall, the results presented 
in Fig. UJ (see also Table [3} allow us to set a rather strong prior 
on the parameter describing the evolution of Yb (see Eq. [4j, with 
—0.02 ^ ay b ^ 0.07, for a conservative range of variation hold- 
ing both at R500 and R2500, accounting for the difference between 
different physical models and for the uncertainties in the estimate of 
the mean associated to the measured intrinsic scatter. As for the nor- 
malization, a conservative allowed range of variation can be taken 
to be Yo,b — 0.83 ± 0.06 at -R2500, with a narrower uncertainty of 
AY ()i i, = 0.03 at i?5oo- 



5 SUMMARY AND DISCUSSION 

In the present study, we have analysed a set of cosmological hydro- 
dynamical simulations of galaxy clusters paying special attention to 
the effects that different implementations of the baryonic physics 
have on the baryon content of these systems. Using th e newest 
versio n of the parallel Tree-PM SPH code GADGET- 3 jSpringell 
120051) . we carried out re-simulations of 29 Lagrangian regions ex- 
tracted around as many galaxy clusters identified within a low- 
resolution N-body parent simulation. These cluster re-simulations 
have been performed using different prescriptions for the baryonic 
physics: without including any radiative processes (NR runs), in- 
cluding the effect of cooling, star formation, SN feedback (CSF 
runs), and including also an additional contribution from AGN 
feedback (AGN runs). 
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Figure 7. The redshift dependence of the mean values of gas depletion Y g (triangles) and baryon depletion Y\, (circles), computed at R500 (left panel) and 
^2500 (right panel), for all clusters that at each redshift have mass M500 > 2 X 10 14 /i — 1 Mq. In each panel, continuous black, dashed-dot blue, and long- 
dashed red lines stand for the NR, CSF, and AGN simulation sets, respectively. These lines have been slightly displaced along the x-axis to avoid overlapping 
among them. Error bars represent la intrinsic scatter computed over all simulated clusters. 



The final sample of objects obtained within each one of these 
set of re-simulations consists in ~ 160 galaxy clusters and groups 
with M vir )3x 10 13 h~ 1 M Q at 2 = 0. Using these three sets 
of simulated galaxy clusters, we have analysed how the different 
physical conditions within them affect their baryon, gas and stellar 
mass fractions and how these results compare with observations. 
We have also examined which are the implications of our results on 
the systematics that affect the constraints on the cosmological pa- 
rameters obtained through the evolution of the cluster baryon mass 
fraction. In order to do so we have calibrated the different bary- 
onic depletion factors and we have analysed their dependences on 
redshift, baryon physics, and cluster radius. 

Our main results can be summarised as follows. 

• In the NR simulations the baryon mass fractions within -R500 
appear flat as a function of the total mass and differ by less than 
~ 10 per cent from the assumed cosmic baryon fraction. This 
result is consis t ent with previous non-radiative simula tions (e.g., 
lEke et al j[l998l ; IKravtsov et alj|2005l : lEttori et alJkoOot) . Whereas 
the CSF simulations present a similar behaviour, when AGN feed- 
back is included there is a significant baryon depletion in poor 
clusters and groups, whereas the cosmic value holds only for the 
most massive clusters. This result, which is in ag reement with th e 
trend displayed by th e obs ervational sam p les from lLin et ail J2003b . 
iGiodini et all | |2009|) . and lLagana e"t~aH teOllh . highlights the effi- 
ciency of the AGN heating in displacing large amounts of gas out- 
side the potential wells of small clusters and groups. 

• The stellar mass fractions obtained in our radiative runs, both 
CSF and AGN, decreases smoothly with increasing cluster mass and 
shows a flattening in the low-mass end « 10 14 M Q ) of our sam- 
ple. When comparing with observational data, the obtained stellar 
mass fractions in our CSF run s is quite large, especia ll y for massive 
cluste rs (e.g. Ilin et al] 120031 : iGonzalez et al.l 120071 ; lLagana et al .1 
l201ll) . When AGN feedback is included, the stellar mass fractions 
within R500 are lowered by about one third, thus alleviating the ten- 
sion with observations, especially at the scale of intermediate-mass 
clusters. However, the level of agreement with observations de- 
pends on the observational sample we compare with. Whereas none 
of our runs is able to reproduce the observed stro ng trend of the 
stellar mass fraction with cluster mass reported bv lGonzalez et al.l 



J2007I) . our results for the AGN runs are in closer agreement with 
other observational samples (e.g jLagana et al.ll201ll) . 

• When analysing the different stellar components separately, 
we find that the fraction of stars within -R500 found in the 
BCG+ICL components is, in the CSF runs, of about 60 per cent 
for massive clusters, and of about 90 per cent for groups. Paradox- 
ically, when AGN feedback is included we find a slightly larger 
fraction of stars in the BCG+ICL component. The reason for this 
is that, although the stellar mass of the BCG+ICL components de- 
creases, the total stellar mass decreases m ore strongly. Thi s resul t, 
in agreement with the AGN simulations bv lPuchwein et al] d2010h . 
is due to the combination of two different effects: while the AGN 
heating truncates the star formation at high redshift, mergers keep 
taking place and unbinding stars from galaxies into the diffuse com- 
ponent, then resulting in a net increase of the fraction of stars that 
end up in the ICL component. 

As for the comparison with observational data, our results con- 
firm a decreasing trend with cluster mass of the fraction of stars 
contributed by BCG and ICL, although they pr edict values of 
this fr action that are larger than those reported bv lGonzalez et al] 
d2007t) , especially for massive clusters. 

• Both of our radiative simulations show that the gas mass 
fraction within clusters increases with increasing cluster mass, 
(see a lso IKravtsov et al ] |2005t iFabian et al.|[2oTol ; Puc hwein et al .1 
l2O10l) . Our results on the values of f g at the scale of rich and 
poor clusters, especially in the presence of AGN feedback, are in 
line with some observational data sets (e.g. , Vikhlini n et al]|200r3 : 
lArnaud et al]| 20071 ; ISun et alj|200°tlZhang et al]|201 lh suggesting 
that low-mass systems have proportionally less gas than high-mass 
systems. 

• The results of our analysis can be used to set priors on the 
parameters which determine the amount of gas depletion when de- 
riving cosmological parameters from observations of the baryon 
fraction in massive clusters. The baryon depletion, Yb, regardless 
of the considered radius or physics, is nearly constant with redshift, 
at least since redshift 2 = 1. However, whereas the obtained evo- 
lution for Yb within R500 is virtually independent of the physics 
included, it shows some dependence on such physical processes 
when looking into inner cluster regions (i?250o). 

• Our results allow us to set a rather strong prior on the pa- 
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rameter describing the evolution of Yj, (see Eq.[4j, with —0.02 ^ 
ctY b ^ 0.07, for a conservative range of variation holding both 
at ifeoo and ifesoo, accounting for the difference between differ- 
ent physical models and for the uncertainties in the estimate of the 
mean associated to the measured intrinsic scatter. As for the nor- 
malization, a conservative allowed range of variation can be taken 
to be Yb,6 = 0.83 ± 0.06 at -R2500, with a narrower uncertainty of 
AY ,i> = 0.03 at R 500 . 

In general, our results show that the star formation in our ra- 
diative runs without AGN heating, even in the presence of rather 
strong galactic winds, is still too efficient, especially in small 
clusters and groups. The situation is significantly improved when 
AGN feedback is included, being able to partly prevent overcool- 
ing in central cluster regions. However, a number of discrepancies 
between simulated and observed baryonic mass fractions within 
clusters still exist, especially when comparing stellar mass frac- 
tions. Nevertheles s, as already reported by other authors (e.g., 
iFabian et alj|20ich . we can infer from our results that a feedback 
source associated to gas accretion onto super-massive BHs seems 
to go in the right direction to conciliate simulations with observa- 
tions. 

Overall, even when the real picture is far more complicated, 
with a number of complex physical processes cooperating to make 
AGN feedback a self-regulated process, we point out that the AGN 
feedback prescription used in the present work significantly im- 
proves previous results on the baryon census within clusters and 
bring closer simulations and observations. 

In general our results highlight that a robust calibration of the 
baryon bias can be defined from simulations at Rsoo> which is quite 
constant within the range of physical models for the ICM included 
in our simulations. This result does not extend at the smaller radius 
^2500, which is the typical radius within which precise measure- 
ments for the gas mass fraction h ave been carried o ut so far for dis- 
tant clusters, using Chandra data j Allen et alj201 ll . and references 
therein). While being beyond the reach of the current generation 
of X-ray telescopes, tracing the gas content of galaxy clusters out 
to large radii requires the next generation of X-ray telescopes to 
be characterized at the same time by a large collecting area and an 
excellent control of the background. 
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Table 1. Best-fit functional forms for the baryon (ft), gas (f g ) and stellar (/*) mass fractions as a function of the total cluster mass, M500, for different 
analyses of observational data. 



Sample 


Best fit 




Lin et al. (2003) 
Giodini et al. (2009) 
Lagana et al. (201 1) 


/b,500 = 
/b,500 = 
/b,500 = 


0A48+°-™l(M 500 /[3 x 10 14 M Q ])(°- 14 8±°-° 4 °> 
(0.123 ± 0.003) (M 500 / [2 x 1O 14 M ])(°- O9±0 - 03 ) 

1O (-O.93O±O.O18) (M5oo/1O 14 M0) (O.136±O.O28) 


Z11+S09 
V06+APP07+S09 


/ g ,500 = lO-C^^O-^CMsoo/llO^MsDC-SOi"-^) 
/ g ,50o(/i/0.7) 3 / 2 = (0.093 ±0.002)(M 50 o/[2 x 10 14 M Q ])(°' 21±0 - 03 ) 


Lin et al. (2003) 
Gonzalez et al. (2007) 
Lagana et al. (201 1) 


/*,500 = 
/*,500 = 
/*,500 = 


0.0164 + °-°°™(Af 5 oo/[3 x 10 14 M Q ])-(°- 2 6±°-09> 

10 (7.57±0.08) M|500 -(0.64±0.13) 
l O C-l.S4±O.lO)( Msoo/[1O 14.5 M0]) (-O.36±O.17) 



Table 2. Values of the gas, stellar and baryonic depletion factors (Y g , Y» and Yj,, respectively) for our set of simulated clusters, for the three different physical 
models (NR, CSF, and AGN), computed at -R2500 an d ^500- I n eacn case, values computed at redshifts z = 0, 0.3, 0.5, 0.8 and 1 are reported. We show 
within brackets the values of the intrinsic scatters computed within the ensamble of simulated clusters. Results are shown for the subset of clusters that, at each 
redshift, are more massive than M 500 = 2 X 10 14 h- 1 M Q . 



Simulation z -R2500 ^500 







Y s 




Y b 


Y g 


Y, 


Yb 


NR 


0.0 


0.80 (0.09) 




0.80 (0.09) 


0.84 (0.04) 




0.84 (0.04) 


NR 


0.3 


0.81 (0.08) 




0.81 (0.08) 


0.85 (0.03) 




0.85 (0.03) 


NR 


0.5 


0.78 (0.09) 




0.78 (0.09) 


0.86 (0.03) 




0.86 (0.03) 


NR 


0.8 


0.84 (0.05) 




0.84 (0.05) 


0.87 (0.03) 




0.87 (0.03) 


NR 


1.0 


0.84 (0.06) 




0.84 (0.06) 


0.86 (0.03) 




0.86 (0.03) 


CSF 


0.0 


0.49 (0.08) 


0.34 (0.07) 


0.85 (0.06) 


0.63 (0.04) 


0.21 (0.04) 


0.85 (0.02) 


CSF 


0.3 


0.48 (0.08) 


0.34 (0.06) 


0.84 (0.05) 


0.63 (0.04) 


0.21 (0.03) 


0.86 (0.03) 


CSF 


0.5 


0.48 (0.06) 


0.32 (0.05) 


0.83 (0.04) 


0.64 (0.03) 


0.20 (0.02) 


0.86 (0.02) 


CSF 


0.8 


0.51 (0.06) 


0.31 (0.05) 


0.86 (0.05) 


0.64 (0.03) 


0.19(0.02) 


0.86 (0.02) 


CSF 


1.0 


0.45 (0.06) 


0.33 (0.06) 


0.82 (0.05) 


0.63 (0.04) 


0.19(0.03) 


0.85 (0.02) 


AGN 


0.0 


0.55 (0.10) 


0.21 (0.03) 


0.77 (0.09) 


0.70 (0.05) 


0.13 (0.02) 


0.85 (0.04) 


AGN 


0.3 


0.54 (0.08) 


0.21 (0.04) 


0.77 (0.07) 


0.70 (0.04) 


0.13 (0.02) 


0.85 (0.04) 


AGN 


0.5 


0.55 (0.07) 


0.22 (0.04) 


0.80 (0.08) 


0.70 (0.03) 


0.13 (0.01) 


0.86 (0.02) 


AGN 


0.8 


0.55 (0.06) 


0.21 (0.02) 


0.79 (0.06) 


0.70 (0.03) 


0.13 (0.01) 


0.85 (0.03) 


AGN 


1.0 


0.51 (0.07) 


0.21 (0.05) 


0.77 (0.09) 


0.68 (0.03) 


0.13 (0.02) 


0.84 (0.04) 



Table 3. Best-fit values of the parameters describing the evolution of the gas and baryonic depletions, according to Eq.|4] For each simulation set (NR, CSF, 
and AGN) and radius of interest (-R250O! ^500> ^200. an d Rvir), we show the normalization (Yo i) ar| d slope (aty ) °f me relation, along with their respective 
standard deviations within brackets, as obtained from the \ 2 minimization procedure. 



Simulation 


Radius 






Y , b 




NR 


Rvir 


0.87 (0.02) 


0.00 (0.04) 


0.87 (0.02) 


0.00 (0.04) 


NR 


R200 


0.86 (0.02) 


0.00 (0.04) 


0.86 (0.02) 


0.00 (0.04) 


IIP. 


R-500 


0.85 (0.03) 


0.02 (0.05) 


0.85 (0.03) 


0.02 (0.05) 


NR 


-R2500 


0.79 (0.07) 


0.07(0.12) 


0.79 (0.07) 


0.07(0.12) 


CSF 


Rvir 


0.70 (0.03) 


-0.03 (0.06) 


0.87 (0.02) 


-0.01 (0.03) 


CSF 


R200 


0.68 (0.03) 


0.00 (0.06) 


0.86 (0.02) 


-0.01 (0.03) 


CSF 


RsOO 


0.63 (0.03) 


0.01 (0.08) 


0.86 (0.02) 


0.00 (0.03) 


CSF 


-R25OO 


0.49 (0.06) 


-0.04(0.18) 


0.85 (0.05) 


-0.02 (0.08) 


AGN 


Rvir 


0.76 (0.03) 


-0.04 (0.05) 


0.87 (0.02) 


-0.01 (0.04) 


AGN 


R200 


0.75 (0.03) 


-0.03 (0.05) 


0.87 (0.03) 


-0.01 (0.04) 


AGN 


R500 


0.71 (0.03) 


-0.03 (0.06) 


0.85 (0.03) 


0.00 (0.05) 


AGN 


-R25OO 


0.55 (0.07) 


-0.04(0.18) 


0.78 (0.07) 


0.01 (0.14) 



